Eggshell_review_project

Compiled under R version 4.2.0 (2022-04-22)

WARNING: edit the working directory to your preferred folder.

This document details all analyses performed in R for the study:
Legendre, L. J., S. Choi, and J. A. Clarke. The use of inconsistent terminology for reptile eggshell traits affects the outcome of evolutionary analyses. Journal of Anatomy 241:641–666.

For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.

Loading packages

library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)

Data and tree

• Tree

tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")

• Data

data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match

## character(0)

rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)

Ancestral state reconstruction: discrete trait (soft/semi-rigid/hard)

• Extract data vectors for each coding

thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)

# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))

# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Using character coding defined in this study (new scoring)

• Perform SIMMAP using phytools – 1000 iterations, using AIC to select best model (code modified from Liam Revell, see here)

x<-thisstudy

aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
    d.aic<-aic-min(aic)
    exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010427592  0.0005213796  0.0005213796
## Semi-rigid  0.0005213796 -0.0010427592  0.0005213796
## Soft        0.0005213796  0.0005213796 -0.0010427592
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010972800  0.0004587151  0.0006385650
## Semi-rigid  0.0004587151 -0.0007221278  0.0002634128
## Soft        0.0006385650  0.0002634128 -0.0009019778
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       0.0000000000  0.000000e+00  0.0000000000
## Semi-rigid 0.0065561168 -1.015329e-02  0.0035971777
## Soft       0.0004585074  4.433192e-05 -0.0005028393
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

logL

##        ER       SYM       ARD 
## -56.60363 -56.31510 -44.90060

AIC<-mapply(aic,logL,c(1,3,6))
AIC

##       ER      SYM      ARD 
## 115.2073 118.6302 101.8012

AIC.W<-aic.w(AIC)
AIC.W

##           ER          SYM          ARD 
## 0.0012254128 0.0002213112 0.9985532760

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
##   1   0 999

o1<-make.simmap(tree,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       0.0000000000  0.000000e+00  0.0000000000
## Semi-rigid 0.0065561168 -1.015329e-02  0.0035971777
## Soft       0.0004585074  4.433192e-05 -0.0005028393
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

o2<-make.simmap(tree,x,model="ER",nsim=nsim[1])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010427592  0.0005213796  0.0005213796
## Semi-rigid  0.0005213796 -0.0010427592  0.0005213796
## Soft        0.0005213796  0.0005213796 -0.0010427592
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

treesstudy<-c(o1,o2)

objstudy<-describe.simmap(treesstudy)

plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 208!

=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.

Using character coding from Norell et al. (2020) – ratio scoring

x<-norell

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard   Semi-rigid         Soft
## Hard       -0.001408430  0.000704215  0.000704215
## Semi-rigid  0.000704215 -0.001408430  0.000704215
## Soft        0.000704215  0.000704215 -0.001408430
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013987686  0.0004876351  0.0009111335
## Semi-rigid  0.0004876351 -0.0010907495  0.0006031144
## Soft        0.0009111335  0.0006031144 -0.0015142479
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.007810285 -0.012947744  0.005137459
## Soft       0.001558250  0.001343177 -0.002901427
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

logL

##        ER       SYM       ARD 
## -70.01788 -69.49006 -56.79592

AIC<-mapply(aic,logL,c(1,3,6))
AIC

##       ER      SYM      ARD 
## 142.0358 144.9801 125.5918

AIC.W<-aic.w(AIC)
AIC.W

##           ER          SYM          ARD 
## 2.685989e-04 6.162366e-05 9.996698e-01

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##   ER  SYM  ARD 
##    0    0 1000

treesnorell<-make.simmap(tree,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.007810285 -0.012947744  0.005137459
## Soft       0.001558250  0.001343177 -0.002901427
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

objnorell<-describe.simmap(treesnorell)

plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Completely different result – let us check which taxa are different from the previous ASR…

datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]

# Only 7 taxa are different – two non-avian dinosaurs and five turtles

# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0014260407  0.0007130204  0.0007130204
## Semi-rigid  0.0007130204 -0.0014260407  0.0007130204
## Soft        0.0007130204  0.0007130204 -0.0014260407
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013980396  0.0005877562  0.0008102834
## Semi-rigid  0.0005877562 -0.0013213236  0.0007335675
## Soft        0.0008102834  0.0007335675 -0.0015438509
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard   Semi-rigid         Soft
## Hard       0.0000000000  0.000000000  0.000000000
## Semi-rigid 0.0076499485 -0.012737081  0.005087133
## Soft       0.0005584858  0.000722445 -0.001280931
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##   ER  SYM  ARD 
##    0    0 1000

treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard   Semi-rigid         Soft
## Hard       0.0000000000  0.000000000  0.000000000
## Semi-rigid 0.0076499485 -0.012737081  0.005087133
## Soft       0.0005584858  0.000722445 -0.001280931
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)

All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.

Using character coding from Legendre et al. (2020) – shell unit scoring

x<-legendre

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007530010  0.0003765005  0.0003765005
## Semi-rigid  0.0003765005 -0.0007530010  0.0003765005
## Soft        0.0003765005  0.0003765005 -0.0007530010
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006161968  0.0000000000  0.0006161968
## Semi-rigid  0.0000000000 -0.0005646245  0.0005646245
## Soft        0.0006161968  0.0005646245 -0.0011808213
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004310504  0.000000000  0.0004310504
## Semi-rigid  0.0036277349 -0.007172395  0.0035446601
## Soft        0.0003893632  0.000000000 -0.0003893632
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

logL

##        ER       SYM       ARD 
## -45.62898 -43.24564 -41.32700

AIC<-mapply(aic,logL,c(1,3,6))
AIC

##       ER      SYM      ARD 
## 93.25795 92.49129 94.65400

AIC.W<-aic.w(AIC)
AIC.W

##        ER       SYM       ARD 
## 0.3372989 0.4948726 0.1678286

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
## 337 495 168

l1<-make.simmap(tree,x,model="ER",nsim=nsim[1])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007530010  0.0003765005  0.0003765005
## Semi-rigid  0.0003765005 -0.0007530010  0.0003765005
## Soft        0.0003765005  0.0003765005 -0.0007530010
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

l2<-make.simmap(tree,x,model="SYM",nsim=nsim[2])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006161968  0.0000000000  0.0006161968
## Semi-rigid  0.0000000000 -0.0005646245  0.0005646245
## Soft        0.0006161968  0.0005646245 -0.0011808213
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

l3<-make.simmap(tree,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004310504  0.000000000  0.0004310504
## Semi-rigid  0.0036277349 -0.007172395  0.0035446601
## Soft        0.0003893632  0.000000000 -0.0003893632
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

treeslegendre<-c(l1,l2,l3)

objlegendre<-describe.simmap(treeslegendre)

plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.

The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.

To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.

• Remove all taxa with differences in coding, except the two sauropodomorphs

remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))

# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Replicate analysis on the new tree 3 times, with identical coding for all taxa except the two sauropodomorphs

• Coded as semi-rigid

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008859112  0.0004429556  0.0004429556
## Semi-rigid  0.0004429556 -0.0008859112  0.0004429556
## Soft        0.0004429556  0.0004429556 -0.0008859112
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -7.746782e-04  4.816389e-05  0.0007265143
## Semi-rigid  4.816389e-05 -7.860575e-04  0.0007378936
## Soft        7.265143e-04  7.378936e-04 -0.0014644079
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001530468  0.0000000000  0.0001530468
## Semi-rigid  0.0042378920 -0.0042378920  0.0000000000
## Soft        0.0020047508  0.0008991591 -0.0029039099
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
## 233 121 646

s1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008859112  0.0004429556  0.0004429556
## Semi-rigid  0.0004429556 -0.0008859112  0.0004429556
## Soft        0.0004429556  0.0004429556 -0.0008859112
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

s2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -7.746782e-04  4.816389e-05  0.0007265143
## Semi-rigid  4.816389e-05 -7.860575e-04  0.0007378936
## Soft        7.265143e-04  7.378936e-04 -0.0014644079
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

s3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001530468  0.0000000000  0.0001530468
## Semi-rigid  0.0042378920 -0.0042378920  0.0000000000
## Soft        0.0020047508  0.0008991591 -0.0029039099
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

treenewsr<-c(s1,s2,s3)

objsr<-describe.simmap(treenewsr)

plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.

• Coded as soft

x[c(21,22)]<-"Soft"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008789375  0.0004394688  0.0004394688
## Semi-rigid  0.0004394688 -0.0008789375  0.0004394688
## Soft        0.0004394688  0.0004394688 -0.0008789375
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       -0.000754388  0.0000000000  0.0007543880
## Semi-rigid  0.000000000 -0.0005361558  0.0005361558
## Soft        0.000754388  0.0005361558 -0.0012905438
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001518191  0.0000000000  0.0001518191
## Semi-rigid  0.0024414178 -0.0024414178  0.0000000000
## Soft        0.0021283836  0.0005467983 -0.0026751820
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
##  59 151 790

so1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008789375  0.0004394688  0.0004394688
## Semi-rigid  0.0004394688 -0.0008789375  0.0004394688
## Soft        0.0004394688  0.0004394688 -0.0008789375
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

so2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       -0.000754388  0.0000000000  0.0007543880
## Semi-rigid  0.000000000 -0.0005361558  0.0005361558
## Soft        0.000754388  0.0005361558 -0.0012905438
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

so3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001518191  0.0000000000  0.0001518191
## Semi-rigid  0.0024414178 -0.0024414178  0.0000000000
## Soft        0.0021283836  0.0005467983 -0.0026751820
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

treenews<-c(so1, so2, so3)

objs<-describe.simmap(treenews)

plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Same result as with semi-rigid coding, unsurprisingly.

• Coded as hard

x[c(21,22)]<-"Hard"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007755266  0.0003877633  0.0003877633
## Semi-rigid  0.0003877633 -0.0007755266  0.0003877633
## Soft        0.0003877633  0.0003877633 -0.0007755266
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006356856  0.000000000  0.0006356856
## Semi-rigid  0.0000000000 -0.000568657  0.0005686570
## Soft        0.0006356856  0.000568657 -0.0012043425
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004481032  0.000000000  0.0004481032
## Semi-rigid  0.0036159226 -0.007150659  0.0035347368
## Soft        0.0003936984  0.000000000 -0.0003936984
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
## 350 494 156

h1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007755266  0.0003877633  0.0003877633
## Semi-rigid  0.0003877633 -0.0007755266  0.0003877633
## Soft        0.0003877633  0.0003877633 -0.0007755266
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

h2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006356856  0.000000000  0.0006356856
## Semi-rigid  0.0000000000 -0.000568657  0.0005686570
## Soft        0.0006356856  0.000568657 -0.0012043425
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

h3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004481032  0.000000000  0.0004481032
## Semi-rigid  0.0036159226 -0.007150659  0.0035347368
## Soft        0.0003936984  0.000000000 -0.0003936984
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333

## Done.

treenewh<-c(h1, h2, h3)

objh<-describe.simmap(treenewh)

plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.

To check how strongly are these results influenced by branch length information, we can replicate these analyses using maximum parsimony, which does not consider branch length information, using castor.

• Coding from this study

studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)

MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

# ancestral states
MPS2<-MPS$ancestral_likelihoods
rownames(MPS2)<-c(209:374); colnames(MPS2)<-c("soft","semi-rigid","hard")

• Coding from Norell et al. (2020)

norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)

MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)

# ancestral states
MPN2<-MPN$ancestral_likelihoods
rownames(MPN2)<-c(209:374); colnames(MPN2)<-c("soft","semi-rigid","hard")

• Coding from Legendre et al. (2020)

legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)

MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

# ancestral states
MPL2<-MPL$ancestral_likelihoods
rownames(MPL2)<-c(209:374); colnames(MPL2)<-c("soft","semi-rigid","hard")

• With the reduced tree, changing only the two sauropodomorphs

xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)

# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

Ancestral state reconstruction – continuous character (eggshell thickness)

• Data and tree

treedi<-multi2di(tree, random=FALSE)

# To visualize node numbers on the tree
ggtree(treedi) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)

• Phylogenetic signal (contMap assumes a Brownian motion model)

phylosig(treedi, AET, method="lambda", test=TRUE)

## 
## Phylogenetic signal lambda : 0.999934 
## logL(lambda) : -315.228 
## LR(lambda=0) : 257.077 
## P-value (based on LR test) : 7.44311e-58

phylosig(treedi, RET, method="lambda", test=TRUE)

## 
## Phylogenetic signal lambda : 0.999934 
## logL(lambda) : -231.158 
## LR(lambda=0) : 207.137 
## P-value (based on LR test) : 5.78775e-47

Very strong signal for both absolute and relative eggshell thickness.

• ASR for absolute eggshell thickness

phyloplotAET<-contMap(treedi, AET, plot=FALSE)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, AET))

## Ancestral character estimates using fastAnc:
##         209         210         211         212         213         214 
##   29.293763   29.293763   26.976449   20.021097   37.072997   57.381022 
##         215         216         217         218         219         220 
##   59.447111   38.957535   37.945828   62.327061   63.691479   60.852431 
##         221         222         223         224         225         226 
##   18.697740   17.464853   16.996346   21.391647   14.489654   17.150500 
##         227         228         229         230         231         232 
##   16.912464   16.433291   14.515531   16.129348    8.218623    8.892389 
##         233         234         235         236         237         238 
##   17.599849   17.653945   11.371660   11.392832   11.093155    7.684835 
##         239         240         241         242         243         244 
##   18.184274   27.793329   11.641005    6.394956    5.262023   11.870736 
##         245         246         247         248         249         250 
##   10.875623   12.261946   12.123157   14.012763   13.895692   14.226702 
##         251         252         253         254         255         256 
##   16.557906   31.172542   39.050721  131.436576  544.304122  573.403854 
##         257         258         259         260         261         262 
##  610.641866  607.826064  498.405442  134.897657  140.865662  142.716126 
##         263         264         265         266         267         268 
##  139.304069  142.860600   61.041297  154.746794  178.179127  198.542856 
##         269         270         271         272         273         274 
##  196.759526  167.261430  145.626598  217.210263  259.744683  324.035230 
##         275         276         277         278         279         280 
##   36.490983  252.266971  278.041830  413.790201  623.886942  366.312141 
##         281         282         283         284         285         286 
##  459.981899  450.032672  468.579008  478.030338   31.529177    2.910210 
##         287         288         289         290         291         292 
##    0.014539    4.051692   11.279474    3.368817   38.203136  115.277972 
##         293         294         295         296         297         298 
##  115.277972 1422.220619 1422.220619  115.277972 1810.898354 1810.898354 
##         299         300         301         302         303         304 
## 1423.663197 1481.546260   29.651406   14.412124   72.463867    1.797803 
##         305         306         307         308         309         310 
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 
##         311         312         313         314         315         316 
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 
##         317         318         319         320         321         322 
## 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 1725.846625 
##         323         324         325         326         327         328 
##  319.220006  479.414045  856.300831  856.300831  856.300831  535.750860 
##         329         330         331         332         333         334 
##  538.224713  538.224713  560.814605 1288.673744  571.363791  610.271423 
##         335         336         337         338         339         340 
##  610.271423  610.271423  610.271423 1392.053831 1101.924581  544.995844 
##         341         342         343         344         345         346 
##  544.995844  544.995844  681.611306  681.611306  420.284152  420.284152 
##         347         348         349         350         351         352 
##  420.284152  401.473506  401.473506  401.473506  401.473506  385.047899 
##         353         354         355         356         357         358 
##  394.698758  483.723330  483.723330  483.723330  483.723330  483.723330 
##         359         360         361         362         363         364 
##  488.479400  483.723330  672.731467  688.660804  743.447690  750.028948 
##         365         366         367         368         369         370 
##  642.663152 1304.111328  319.380290  778.921709  908.582073  820.930413 
##         371         372         373         374         375         376 
##  423.106033  255.141714  300.722998  343.668052  298.875220  290.389862 
##         377         378         379         380         381         382 
##  222.173714  217.770862  177.872861  166.751051  161.039115  229.707740 
##         383         384         385         386         387         388 
##  234.260007  236.237640  233.565993  232.205729  231.164466  196.488111 
##         389         390         391         392         393         394 
##  216.864203  210.923903  208.771974  210.728366  206.860831  201.759935 
##         395         396         397         398         399         400 
##  251.229188  289.706525  316.999616  253.504459  257.542381  243.366955 
##         401         402         403         404         405         406 
##  151.843383  155.880089  164.029302  144.791426  128.410742  106.962875 
##         407         408         409         410         411         412 
##  100.299131  118.844349  109.671081   73.950228  103.670078   99.287810 
##         413         414         415 
##   98.491772  103.798544  101.498594

Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.

Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).

• ASR for relative eggshell thickness

phyloplotRET<-contMap(treedi, RET, plot=FALSE)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, RET))

## Ancestral character estimates using fastAnc:
##       209       210       211       212       213       214       215       216 
##  3.435293  3.435293  4.047709  7.765196 20.362606 37.938372 39.824787 25.690342 
##       217       218       219       220       221       222       223       224 
## 39.657428 42.678094 44.856834 53.006727  7.508812  7.363644  7.173682  6.553561 
##       225       226       227       228       229       230       231       232 
##  9.463168  7.383174  7.655193  9.189493 11.648581 13.022670 12.734649 13.646320 
##       233       234       235       236       237       238       239       240 
## 16.530771 20.145434 21.874698 24.370253 23.820866 14.769114  1.647011  1.003134 
##       241       242       243       244       245       246       247       248 
##  1.610621  0.936574  0.860310  1.828622  1.596660  2.012543  3.474795  1.682980 
##       249       250       251       252       253       254       255       256 
##  1.180146  1.164828  2.085553  3.133934  2.630136  6.622457 10.830493 10.810476 
##       257       258       259       260       261       262       263       264 
## 10.535636 10.641238  9.256621  6.895434  8.089943 10.710668 10.792255 12.171031 
##       265       266       267       268       269       270       271       272 
##  2.043677 11.100157 10.747658 10.935592 10.272658  9.993139  8.767617 10.221160 
##       273       274       275       276       277       278       279       280 
##  7.286643  7.960715  2.157566  4.753714  4.814619  5.477856  6.779844  4.588011 
##       281       282       283       284       285       286       287       288 
##  5.809998  6.065841  5.670088  5.485997  1.720444  0.460659  0.004019  0.552030 
##       289       290       291       292       293       294       295       296 
##  1.153563  0.408696  1.278576  1.232461  1.232461  5.333110  5.333110  1.232461 
##       297       298       299       300       301       302       303       304 
##  0.657841  0.657841  1.223907  1.273789  0.903462  0.431246  0.388596  0.130004 
##       305       306       307       308       309       310       311       312 
##  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681 
##       313       314       315       316       317       318       319       320 
##  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681  0.730681 
##       321       322       323       324       325       326       327       328 
##  0.730681  0.730681  2.766052  4.215664  4.204283  4.204283  4.204283  5.022112 
##       329       330       331       332       333       334       335       336 
##  5.237510  5.237510  5.164118  2.788452  5.172184  4.838514  4.838514  4.838514 
##       337       338       339       340       341       342       343       344 
##  4.838514  2.039981  3.278505  5.533461  5.533461  5.533461  4.218526  4.218526 
##       345       346       347       348       349       350       351       352 
##  8.313752  8.313752  8.313752  8.642456  8.642456  8.642456  8.642456  8.868480 
##       353       354       355       356       357       358       359       360 
##  8.729123  8.015217  8.015217  8.015217  8.015217  8.015217  8.060607  8.015217 
##       361       362       363       364       365       366       367       368 
##  3.120847  3.014030  2.340587  2.257551  2.382146  0.661596  4.633715  2.110811 
##       369       370       371       372       373       374       375       376 
##  1.859324  1.877740  1.176540 10.823588  9.306410  7.854842  9.850179 10.769445 
##       377       378       379       380       381       382       383       384 
## 11.854208 11.578107 13.301104 13.850234 14.334685 10.815390 10.465704 10.017208 
##       385       386       387       388       389       390       391       392 
##  9.675993  9.353371  9.030021 13.460333 11.098649 10.827129 10.900202 10.571104 
##       393       394       395       396       397       398       399       400 
## 10.467844 10.384180  9.573346  8.310444  7.454699  9.414721  8.811303  9.272529 
##       401       402       403       404       405       406       407       408 
## 19.513793 20.272830 19.450144 20.560878 23.253747 28.368872 30.395795 25.079676 
##       409       410       411       412       413       414       415 
## 27.481017 41.577478 29.372089 30.523206 32.152446 31.107990 34.365793

Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.

We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).

References

• Legendre, L.J., Clarke, J.A., 2021. Shifts in eggshell thickness are related to changes in locomotor ecology in dinosaurs. Evolution 75, 1415–1430.
• Legendre, L.J., Rubilar-Rogers, D., Musser, G.M., Davis, S.N., Otero, R.A., Vargas, A.O., Clarke, J.A., 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.
• Norell, M.A., Wiemann, J., Fabbri, M., Yu, C., Marsicano, C.A., Moore-Nall, A., Varricchio, D.J., Pol, D., Zelenitsky, D.K., 2020. The first dinosaur egg was soft. Nature 583, 406–410.
• Stein, K., Prondvai, E., Huang, T., Baele, J.-M., Sander, P.M., Reisz, R., 2019. Structure and evolutionary implications of the earliest (Sinemurian, Early Jurassic) dinosaur eggs and eggshells. Scientific Reports 9, 4424.